Analizie poddany został zbiór danych dotyczących połowu śledzia atlantyckiego oraz warunków środowiska w jakim żyje. Celem analizy była próba odpowiedzi na pytanie dlaczego od pewnego czasu złowione ryby mają coraz mniejszy rozmiar. Eksploracja rozpoczyna się od załadowania odpowiednich bibliotek, następnie ustawione zostaje ziarno zapewniające powtarzalność uzyskanych wyników oraz wczytany zostaje zbiór danych. Sprawdzona zostaje liczba pustych wartości oraz kształt zbioru. Ze względu na relatywnie niską ilość poszczególnych wartości NA zostaje podjęta decyzja o ich usunięciu.
Później przyglądano się bliżej poszczególnym atrybutom, ich rozkładom wartości oraz postawowym statyskom (mean, min, max, sd, etc.).
Kolejnym krokiem było sprawdzenie wzajemnych zależności pomiędzy atrybutami (włączając w to zmienną przewidywaną). Większość atrybutów nie wykazuje znaczącej wzajemnej zależności. Jednak kilka zmiennych, w szczególności dotyczących dostępności planktonu, jest silnie skorelowane. Dla zmiennych wykazujących wzajemne powiązania dokonano dodatkowej wizualizacji. Żaden z atrybutów nie jest bezpośrednio skorelowany z wartością zmiennej przewidywanej.
Następnie sporządzono wykres zmiany długości śledzia w czasie, jasno wskazujący na spadek średniego rozmiaru złowionych ryb.
Następny etap opierał się na zbudowaniu regresora, mającego za cel predykcję długości śledzia. Dokonano podziału zbioru na część treningową, walidacyjną oraz testową. W pierwszej kolejności przetestowano prosty model regresji liniowej i zmierzono jego jakość predykcji przy pomocy metryk RMSE oraz R2. Drugim wybranym modelem był XGBoost - wydajna biblioteka udostępniająca rozwiązania z dziedziny gradient boosting, oparte na drzewach decyzyjnych. Jest często używany w problemach analizy danych. Model wymagał dodatkowej konwersji danych na kompatybilny format. Następnie przy pomocy zbioru walidacyjnego dobrano parametry, tak, aby zmniejszyć overfitting. Użycie XGBoost pozwoliło na wyraźną poprawę wyników predykcji.
Finalnie, dla modelu o najlepszej skuteczności(XGBoost), przy pomocy metody LIME dokonano analizy ważności atrybutów oraz ich wpływu na estymowaną wartość w poszczególnych przypadkach.
Do najważniejszych cech należą odwiednio:
sst - 0.57,recr - 0.14,lcop1 - 0.06.XGBoost dokonuje automatycznej analizy zależności zmiennych, na wykresie feature importance można zauważyć, że tylko jedna z silnie skorelowanych zmiennych ma wyższe znaczenie dla modelu.
Z dużego wpływu temperatury przy powierzchni wody (sst) na podejmowaną przez model decyzję można by wnioskować, że drobny wzrost temperatury (mniej więcej od podobnego okresu, gdy długość śledzia zaczęła spadać) mógł mieć wpływ na zmniejszenie rozmiaru poławianych śledzi.
library(reshape2)
library(plyr)
library(dplyr)
library(caret)
library(corrplot)
library(ggplot2)
library(plotly)
library(glmnet)
library(gbm)
library(xgboost)
library(lime)
library(nnet)
set.seed(10)
Wartości puste oznaczone są znakiem “?”, przy wczywaniu zbioru korzystamy z parametru na.strings do zastąpienia ich wartością NA
data <- read.csv("sledzie.csv", header=TRUE, na.strings = "?")
nmissing <- function(x) sum(is.na(x))
knitr::kable(colwise(nmissing)(data))
| X | length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | recr | cumf | totaln | sst | sal | xmonth | nao |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 1581 | 1536 | 1555 | 1556 | 1653 | 1591 | 0 | 0 | 0 | 0 | 1584 | 0 | 0 | 0 |
Wymiary zbioru: (52582, 16)
Całkowita liczba wierszy: 52582
Liczba wierszy zawierających pustą wartość: 10094
Liczba wierszy po usunięciu pustych wartości: 42488
| X | length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | recr | cumf | totaln | sst | sal | xmonth | nao | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 24923 | 24922 | 24.0 | 0.00000 | 0.70118 | 11.50000 | 5.67765 | 23.00000 | 9.17171 | 0.227 | 421391 | 0.1480941 | 730351.2 | 13.63160 | 35.50835 | 2 | -2.86 |
| 18378 | 18377 | 25.5 | 0.00000 | 0.72727 | 2.59316 | 33.12121 | 5.35088 | 37.23232 | 0.547 | 148045 | 0.2714717 | 179955.4 | 13.39520 | 35.49992 | 5 | -2.14 |
| 22182 | 22181 | 25.5 | 0.00000 | 0.70118 | 11.50000 | 5.67765 | 23.00000 | 9.17171 | 0.227 | 421391 | 0.1480941 | 730351.2 | 13.63160 | 35.50835 | 4 | -2.86 |
| 34827 | 34826 | 25.0 | 3.14462 | 5.80145 | 4.26148 | 18.43839 | 7.67148 | 27.47359 | 0.254 | 474983 | 0.2073709 | 534157.2 | 13.89253 | 35.50455 | 5 | -0.75 |
| 46082 | 46081 | 22.5 | 0.05455 | 0.23134 | 1.48629 | 11.67093 | 1.83515 | 17.80753 | 0.591 | 473462 | 0.3758273 | 306067.6 | 14.35600 | 35.52329 | 9 | 0.72 |
W zbiorze znajdują się następujące kolumny(atrybuty):
Na histogramie długości złowionego śledzia wyraźnie rysuje się rozkład normalny. Odchylenie standardowe jest relatywnie niskie, równe jest 1,65 cm.
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 19.00 24.00 25.50 25.30 26.50 32.50 1.65
Dostępność planktonu cfin1 w okresie pomiarów była zdecydowanie niska, rozkład wartości jest w ponad 80% zdominowany przez wartość minimalną. Zdarzają się nieliczne wysokie wartości.
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 0.0000 0.0000 0.1111 0.4457 0.3333 37.6667 0.9800
Podobnie jak w przypadku cfin1, rozkład cfin2 jest zdominowany przez niskie wartości, jednak w znacznie mniejszym stopniu. W tym przypadku również występują wyższe wartości cechy, jednak tutaj układają się w widoczne grupy.
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 0.0000 0.2778 0.7012 2.0269 1.7936 19.3958 3.7200
Analogicznie jak pozostałych przypadkach widoczna jest niska dostępność planktonu chel1.
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 0.000 2.469 5.750 10.016 11.500 75.000 14.320
Dostępność chel2 w całym zbiorze jest dość zróżnicowana, w rozkładzie występuje dużo różnych wartości o podobnej częstości występowania.
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 5.238 13.427 21.435 21.197 27.193 57.706 9.990
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 0.3074 2.5479 7.0000 12.8386 21.2315 115.5833 15.1000
Podobnie jak chel2, występowanie planktonu lcop2 jest dość zróżnicowane. Wartości obu planktonów wydają się być skorelowane.
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 7.849 17.808 24.859 28.396 37.232 68.736 13.870
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 0.0680 0.2270 0.3320 0.3306 0.4650 0.8490 0.1600
Zmienna recr przedstawia roczny narybek reprezentowany w postaci liczby złowionych ryb, wartości liczbowe osiągają dość wysokie wartości, niejednokrotnie przekraczające milion.
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 140515 360061 421391 519877 724151 1565890 270639
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 0.06833 0.14809 0.23191 0.22987 0.29803 0.39801 0.09000
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 144137 306068 539558 515082 730351 1015595 221389
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 12.77 13.60 13.86 13.87 14.16 14.73 0.44
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 35.40 35.51 35.51 35.51 35.52 35.61 0.04
Dane dotyczące miesiąca mają w przybliżeniu rozkład normalny, większość pomiarów pochodzi z okresu letnio-jesiennego.
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## 1.000 5.000 8.000 7.252 9.000 12.000 2.760
## Min. 1st Qu. Median Mean 3rd Qu. Max. Std.Dev.
## -4.89000 -1.90000 0.20000 -0.09642 1.63000 5.08000 2.25000
Z wykresu można wywnioskować, że wiekszość zmiennych jest skorelowana ze sobą jedynie w lekkim stopniu lub wcale. Jest jednak kilka wyjątków, przede wszystkim w przypadku dostępności niektórych gatunków planktonu. Występują dwie bardzo silne korelacje:
0.96 między chel1, a lcop1,0.89 między chel2, a lcop2.oraz kolejna para silnych korelacji dodatnich:
0.82 między fbar, a cumf0.65 między cfin2, a lcop2i jedna silna korelacja ujemna:
-0.71 między totaln, a cumfNa wykresie widzimy bardzo silnie skorelowane zmienne lcop1 oraz chel1
W tym przypadku mamy do czynienia z negatywną korelacją zmiennych totaln i cumf. Siła korelacji jest mniejsza niż w poprzednich przypadkach, jednak wciąż można łatwo dostrzec malejący trend.
Zaprezentowany wykres, celem zwiększenia czytelności, jest przygotowany na bazie próbki 1000 pomiarów. Dzięki wygładzonej linii średniej można łatwo dostrzec jaki był kształt trendu w czasie. Interaktywny wykres pozwala dokładnie prześledzić zmianę wartości, na wykresie wyraźnie widać, że po około 16000. próbce średnia długość śledzia zaczyna spadać, dzieje się tak aż do końca badanego okresu.
train_p <- 0.7
valid_p <- 0.5
idx_train <- createDataPartition(y = data_clean$length, p = train_p, list = FALSE)
train_set <- data_clean[,-1][ idx_train,]
test_set <- data_clean[,-1][-idx_train,]
idx_valid <- createDataPartition(y = test_set$length, p = valid_p, list = FALSE)
valid_set <- test_set[ idx_valid,]
test_set <- test_set[-idx_valid,]
Rozmiary zbiorów: treningowego, walidacyjnego oraz testowego
dim(train_set)
## [1] 29744 15
dim(valid_set)
## [1] 6372 15
dim(test_set)
## [1] 6372 15
simple_regr <- lm(length~. ,data = train_set)
simple_prediction <- predict(simple_regr, test_set[,-1])
Podsumowanie parametrów modelu regresji lm
summary(simple_regr)
##
## Call:
## lm(formula = length ~ ., data = train_set)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.8848 -0.9335 0.0018 0.9126 6.9429
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.487e+01 8.458e+00 7.669 1.79e-14 ***
## cfin1 1.384e-01 9.372e-03 14.769 < 2e-16 ***
## cfin2 1.698e-02 5.224e-03 3.250 0.00115 **
## chel1 -1.023e-02 2.250e-03 -4.546 5.50e-06 ***
## chel2 -4.628e-03 3.316e-03 -1.396 0.16284
## lcop1 2.181e-02 2.097e-03 10.401 < 2e-16 ***
## lcop2 1.129e-02 2.756e-03 4.094 4.25e-05 ***
## fbar 6.126e+00 1.132e-01 54.135 < 2e-16 ***
## recr -3.905e-07 3.729e-08 -10.470 < 2e-16 ***
## cumf -1.021e+01 2.112e-01 -48.323 < 2e-16 ***
## totaln -5.909e-07 7.139e-08 -8.277 < 2e-16 ***
## sst -1.274e+00 2.674e-02 -47.659 < 2e-16 ***
## sal -6.086e-01 2.395e-01 -2.542 0.01104 *
## xmonth 9.068e-03 2.863e-03 3.167 0.00154 **
## nao 3.519e-02 5.582e-03 6.305 2.92e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.358 on 29729 degrees of freedom
## Multiple R-squared: 0.3224, Adjusted R-squared: 0.3221
## F-statistic: 1010 on 14 and 29729 DF, p-value: < 2.2e-16
Wyniki
R2 0.3156754
RMSE 1.3586153
Przygotowanie zbioru do wersji kompatybilnej z XGBoost
xgb_train_set <- xgb.DMatrix(data = as.matrix(train_set[,-1]), label = train_set$length)
xgb_valid_set <- xgb.DMatrix(data = as.matrix(valid_set[,-1]), label = valid_set$length)
xgb_test_set <- xgb.DMatrix(data = as.matrix(test_set[,-1]), label = test_set$length)
watchlist <- list(train=xgb_train_set, valid=xgb_valid_set)
xg <- xgb.train(data = xgb_train_set, max.depth=32, nthread = 5, nrounds = 20, watchlist=watchlist)
## [1] train-rmse:17.424788 valid-rmse:17.426558
## [2] train-rmse:12.230176 valid-rmse:12.232804
## [3] train-rmse:8.605362 valid-rmse:8.609152
## [4] train-rmse:6.084421 valid-rmse:6.088549
## [5] train-rmse:4.341674 valid-rmse:4.347968
## [6] train-rmse:3.150393 valid-rmse:3.159253
## [7] train-rmse:2.353688 valid-rmse:2.366677
## [8] train-rmse:1.837287 valid-rmse:1.854480
## [9] train-rmse:1.517758 valid-rmse:1.538785
## [10] train-rmse:1.331451 valid-rmse:1.356127
## [11] train-rmse:1.228172 valid-rmse:1.255594
## [12] train-rmse:1.173326 valid-rmse:1.203040
## [13] train-rmse:1.144919 valid-rmse:1.176587
## [14] train-rmse:1.130379 valid-rmse:1.163402
## [15] train-rmse:1.122979 valid-rmse:1.156788
## [16] train-rmse:1.119246 valid-rmse:1.153730
## [17] train-rmse:1.117372 valid-rmse:1.152286
## [18] train-rmse:1.116414 valid-rmse:1.151726
## [19] train-rmse:1.115926 valid-rmse:1.151528
## [20] train-rmse:1.115675 valid-rmse:1.151489
Wyniki
R2 0.5224293
RMSE 1.1351347
Przy pomocy pakietu Lime można zaobserwować wpływ poszczególnych atrybutów na wartość przewidywanej zmiennej dla kilku wybranych próbek. Analizując kilka wybranych próbek można łatwo zauważyć pewne prawidłowości:
sst > 14.2 bardzo negatywnie wpływa rozmiar śledzia,lcop1 > 21.23 ma lekki pozytywny wpływ,W przypadku większej ilości przypadków można je zwizualizować na wykresie w postaci heatmap.
XGBoost oferuje wbudowane feature importance mierzone przy pomocy parametru gain. Pośród użytych atrybutów zdecydowanie dominuje sst, mając ponad połowę z całego zysku przy podejmowaniu decyzji. Następny jest atrybut recr mający ponaddwukrotnie większą ważność od swojego następnika. Pozostałe zmienne mają już mniejszy, ale bardziej równomierny wkład w decyzję.
importanceRaw <- xgb.importance(feature_names = colnames(test_set[,-1]), model = xg, data = xgb_test_set)
xgb.plot.importance(importance_matrix = importanceRaw)
Z dużego wpływu temperatury przy powierzchni wody (sst) na podejmowaną przez model decyzję można by wnioskować, że drobny wzrost temperatury (mniej więcej od podobnego okresu, gdy długość śledzia zaczęła spadać) mógł mieć wpływ na zmniejszenie rozmiaru poławianych śledzi.